Draft version June 16, 2009 

Preprint typeset using LAT^^ style cmulateapj v. 11/12/01 



SIMULATIONS OF ELECTRON MAGNETOHYDRODYNAMIC TURBULENCE 

JUNGYEON ChO^'^ 
AND 

A. Lazarian^ 
Draft version June 16, 2009 

ABSTRACT 

We present numerical simulations of electron magnetohydrodynamic (EMHD) and electron reduced 
MHD (ERMHD) turbulence. Comparing scaling relations, we find that both EMHD and ERMHD 
turbulence show similar spectra and anisotropy. We develop new techniques to study anisotropy of 
EMHD turbulence. Our detailed study of anisotropy of EMHD turbulence supports our earlier result 

1/3 

of A:|| cx scaling, where fcy and k± are wavenumbers parallel and perpendicular to local direction of 
magnetic field, respectively. We find that the high-order statistics show a scaling that is similar to the She- 
Leveque scaling for incompressible hydrodynamic turbulence and different from that of incompressible 
MHD turbulence. We observe that the bispectra, which characterize the interaction of different scales 
within the turbulence cascade, are very different for EMHD and MHD turbulence. We show that 
both decaying and driven EMHD turbulence have the same statistical properties. We calculate the 
probability distribution functions (PDFs) of MHD and EMHD turbulence and compare them with those 
of interplanetary turbulence. We find that, as in the case of the solar wind, the PDFs of the increments 
of magnetic field strength in MHD and EMHD turbulence are well described by the Tsallis distribution. 
We discuss implications of our results for astrophysical situations, including the advection dominated 
accretion flows and magnetic reconnection. 

Subject headings: MHD — turbulence — solar wind 



1. INTRODUCTION 

Turbulence at scales below the proton gyroradius is of 
great importance in many astrophysical applications. Such 
turbulence involving motions of electrons is essential for 
understanding the small-scale magnetic field dynamics of 
plasmas. It is also important for understanding of mag- 
netic fields in the crust of a neutron star (Goldreich & 
Reisenegger 1992). This turbulence has been measured at 
by solar wind probes (Leamon et al. 1998). 

The origin of the small-scale turbulence in a magnetized 
plasma is easy to understand if we think of what is happen- 
ing with turbulent motions at small scales. When turbu- 
lence is driven on large scales, turbulence energy cascades 
down to smaller scales. The nature of magnetized turbu- 
lence from the outer scale to the proton gyroradius scale 
is relatively well known. Magnetized turbulence above 
the proton gyroradius can be described by the standard 
magnetohydrodynamics (MHD). As the name implies, the 
standard MHD treats the plasma as a single fluid. MHD 
turbulence can be decomposed into cascades of Alfven, 
fast and slow modes (Goldreich & Sridhar 1995; Lithwick 
& Goldreich 2001; Cho & Lazarian 2002, 2003). While 
fast and slow modes get damped at larger scales, Alfvenic 
modes can cascade down to the proton gyroradius scale. 
Near and below the proton gyroradius scale, single-fluid 
description fails and we should take into account kinetic 
effects. Then what will happen to the Alfven modes that 
reach the proton gyroradius scale? 

Recent years, the nature of small-scale MHD turbu- 
lence in the solar wind has drawn interest from researchers 
(Howes et al. 2008a, b; Matthaeus, Servidio, & Dmitruk 



2008; Saito et al. 2008; Gary, Saito, & Li 2008; Schekochi- 
hin et al. 2009; Dmitruk & Matthaeus 2006). In situ mea- 
surements of the solar wind show magnetic fluctuations 
over a broad range of frequencies. In the rest frame of 
the spacecrafts, the magnetic fluctuations show a broken 
power-law spectrum. For example, Leamon et al. (1999) 
reported that, at ^ 0.2Hz, the spectrum breaks from a 
power-law to a i^^^-^^ power-law. In general, the 
spectral break-point lies in the range 0.2Hz ^v< 0.5Hz 
(see Saito et al. 2008). The i^^^-^'^ power-law ati^ < 0.2Hz 
seems to be relatively robust and represents inertial range 
of Alfvenic MHD turbulence. However, the power index 
for v > 0.5Hz seems to vary between -2 and -4 (Leamon 
et al. 1998; Smith et al. 2006). This range, character- 
ized by a steeper power-law index, is termed "dispersion 
range" (Stawicki, Gary, & Li 2001), which is different from 
the dissipation range. The true dissipation scale of turbu- 
lence may lie at the end of the dispersion range. When 
we convert frequency to length scale, the broken power- 
law implies that the magnetic energy spectrum changes 
from a inertial range spectrum to a steeper disper- 

sion range spectrum, as we move from large scales to small 
scales. The transition from the inertial range to the dis- 
persion range occurs near the proton gyro-scale pi . (How- 
ever, it is also possible that it occurs at the ion inertial 
scale di = Pi/^/^, where (3i is the ion plasma /3. See dis- 
cussion in Schekochihin et al. 2009). The identity of the 
dispersion range turbulence is still under debate. 

Small-scale magnetized turbulence also plays important 
roles in other astrophysical objects. Among them, the 
crust of neutron stars gives us useful insights on the elec- 
tron MHD (EMHD) model of small-scale magnetized tur- 



^ Dept. of Astronomy and Space Science, Chungnam National Univ., Daejeon, Korea; cho@canopus.cnu.ac.kr 
^ Dept. of Astronomy, Univ. of Wisconsin, Madison, WI53706, USA; Iazarian@astro.wisc.edu 



2 



Cho & Lazarian 



bulence. In the crust of neutron stars, ions are virtu- 
ally immobile and, thus, the ion gyroradius scale can be 
regarded as infinite. Therefore, turbulence in the crust 
should be similar to small-scale turbulence. Since ions are 
immobile, they provide a smooth charge background and 
electrons carry all the current, so that 



-V X B, 



where Ve is the electron velocity, J is electric current den- 
sity, B is magnetic field, c is the speed of light, rig is the 
electron number density, and e is the absolute value of the 
electric charge. Inserting this relation into the magnetic 
induction equation, we can obtain the EMHD equation 

5B 



dt 



ATTUpe 



V X [(V X B) X B] ryV^B, (2) 



where rj is magnetic diffusivity (see Kingscp, Chukbar, 
& Yankov 1990 for details about EMHD). Goldreich & 
Reisenegger (1992) first showed that magnetized turbu- 
lence in the crust of neutron stars can be described by the 
EMHD equation. They discussed the properties of EMHD 
turbulence in neutron stars and argued that EMHD tur- 
bulence can enhance ohmic dissipation of magnetic field 
in isolated neutron stars (see also Gumming, Arras, & 
Zweibel 2004). 

In this paper, we will focus on spectrum and anisotropy 
of EMHD turbulence. Earlier researchers convincingly 
showed that energy spectrum of EMHD turbulence is 
steeper than Kolmogorov's k~^/^ spectrum (Biskamp, 
Schwarz, & Drake 1996; Biskamp et al. 1999; Ng et 
al.2003). They found that the energy spectrum follows 



E{k) oc fc"^/^ 



(3) 



turbulence. In §4, we describe our numerical methods. In 
§5, we compare EMHD and ERMHD turbulence. In the 
section we perform simulations in an elongated numeri- 
cal box (786 X 256^). In §6, we present detailed study of 
anisotropy from an EMHD simulation with 512'^ resolu- 
tion. In §7, we present high-order statistics and bispectra 
of EMHD turbulence. In §8, we calculate the probability 
distribution functions (PDFs) and briefly compare them 
with the solar wind data. We give discussions in §9 and 
summary in §10. 

2. EMHD AND ERMHD 

Alfven modes are incompressible and thus less prone 
to collisionless damping (Barnes 1966; Kulsrud & Pearse 
1969). Therefore, in the solar wind, it is generally accepted 
that Alfvenic turbulence provides a suitable description of 
fluid motions on scales larger than the proton gyroscale. 
The physics of strong Alfven turbulence is relatively well 
understood (see, e.g., Goldreich & Sridhar 1995) 

However, turbulence on scales smaller than the proton 
gyroscale is not well understood. Nevertheless, there are 
several models for this small-scale turbulence. Here, we 
consider only fluid-like models of the small-scale turbu- 
lence. In particular, we consider Electron MHD (EMHD) 
model and Electron Reduced MHD (ERMHD) model. The 
former has been studied since 1990s (Kingsep et al. 1990; 
Biskamp et al. 1996) and the latter has been proposed 
only recently (Schekochihin et al. 2009). Both models can 
be derived from the generalized Ohm's law. 

2.1. EMHD 

EMHD can be viewed HaU MHD in the hmit of fcp^ > 1, 
where pi is the ion gyroradius. Hall MHD equations are 
very similar to the standard MHD ones. But, there are 
also differences. The most important difference is that 
Hall MHD is based on the generalized Ohm's law 

J X B J 



E = X B 

c 



The steep energy spectrum can be explained by the fol- 
lowing Kolmogorov-type argument (Biskamp et al. 1996). 
Suppose that the eddy interaction time for eddies of size 
I is the usual eddy turnover time teas,/ ^ l/vi- Since 
V oc V X B (Eq. [1]), this becomes teas,/ oc P/bi. Gombin- 
ing this with the constancy of spectral energy cascade rate 
{bf /tcas,i=constant), one obtains E{k) oc k~'^/^. Note that 
E(k) and hi are related by kE{k) ~ 6^. Since earlier re- 
searchers convincingly obtained the k~'^/'^ spectrum, more 
focus is given to anisotropy of EMHD turbulence. 

Gho & Lazarian (2004; hereinafter GL04) derived the 
expression for anisotropy 

/cj| oc k^^ (4) 

where fcj^ and /c|| should be understood as wavenumbers 
parallel and perpendicular to the local magnetic field, the 
same way as parallel and perpendicular wavenumbers are 
understood in Goldreich- Sridhar (1995) model of MHD 
turbulenc^. The expression (|4]), as we discuss in §3, fol- 
lows from the application of the critical balance notion to 
the electron MHD cascade. 

In this paper, we only consider strong EMHD tur- 
bulence. Discussions on weak EMHD turbulence can 
be found in Galtier & Bhattacharjee (2003) and Galtier 
(2006). In §2, we compare EMHD and recently proposed 
ERMHD (Schekochihin et al. 2009) formalisms. In §3, we 
discuss expected scaling relations of EMHD and ERMHD 

A wavelet description would be more appropriate, as usual wavenumbers are defined in the global magnetic field frame. We keep this in mind, 
while using the traditional wavenumber notations. 



(5) 

while the standard MHD uses the 'standard' Ohm's law 

V J 

E = xB + -. (6) 

C cr 

Here, v is fluid velocity, E is electric field, and a is con- 
ductivity. Gompared with that of the standard MHD, the 
induction equation of Hall MHD has an extra term (the 
J X B term on the right-hand side): 
SB „ ^ 



V X (v X B) - V X 

V X (v X B) - V X 



J X B 



c(V X B) X B 



-V^B 



TyV^B, (7) 



47r rige 

where -q = /{Ana) and we use V x B = (47r/c)J. 

When we normalize magnetic field to a velocity 
(i.e. B/^47rp 



/ ^ AT^niirii — s- B, where rii is the ion 
number density and is the ion mass), Eq. ([7]) becomes 
(}B 

— = Vx(vxB)-djVx[(V X B) X B]-H?7V^B, (HaU MHD) 

(8) 



EMHD Turbulence 



3 



where 



■iTTUee 



Pi 



(9) 



yj'4:Trni{Ze)^ /rrii ^^pi \fWi 

Here, pi is the ion gyroradius, Z = n(./ni = Qi/e, and /3, 
is the ion plasma beta: 

UikBT 



|3^ 



(10) 



= -d,V X [(V X B) X B] + T^V^B. (EMHD) (13) 



B2/87r' 

where ks is the Boltzmann constant. The order of magni- 
tude values of the first and the second terms on the right- 
hand side are 

V X (v X B) - h'^/l, (11) 

(i,V X [(V X B) X B] - d^h^/f, (12) 

where I is the scale of interest and we assume v ^ b. 
Eq. ([5]) reduces to the standard MHD induction equation 
for I ^ di, while it reduces to the Electron MHD equation 
for I <^ di'. 

dB 

'dt 

In usual coUisionless plasmas, transition from standard 
MHD to EMHD occurs at the ion incrtial scale di. When 
ions arc immobile and provide a homogeneous background, 
as in the crust of a neutron star, EMHD can be applicable 
even at the outer scale of turbulence, which can be larger 
than di defined by the first two expressions in Eq. 
Note that anisotropy of turbulence (i.e. k± 3> fc||) at di 
(in usual plasmas) or on the energy injection scale (in case 
ions are immobile) is not a necessary condition for EMHD. 

2.2. ERMHD 

Schekochihin et al. (2009) first derived ERMHD from ki- 
netic RMHD equations. They also gave derivation of ER- 
MHD equations from the generalized Ohm's law. There- 
fore the starting point of ERMHD may be also the gener- 
alized Ohm's law and derivation of ERMHD is identical to 
the EMHD case up to Eq. ([7]) (sec previous subsection). 
Note that the velocity v in Eq. ([7]) denotes ion velocity 
and that RMHD, hence ERMHD, assumes anisotropy of 
turbulence (i.e. fc_L S> fc||). 

However, ERMHD assumes that the term 

V X (v X B) « -BV • V (14) 

in Eq. ([7]) may not be negligible in the limit of kpi 3> 1. 
That is, ERMHD assumes w in this limit, but 
V • Vi = V • Ve 7^ (see Schekochihin et al. 2009). Taking 
this into account, one can rewrite the magnetic induction 
equation as 

SB _ ^ c_ 

dt * Airen, 

which becomes 



■V X [B • VB], 



(15) 



SB, I 



9B^ _ 
-BoV • V, 



47ren, 



■V_L X [B • VB||], (16) 
-Vix[B-VBJ, (17) 



dt " ' AneUe 

where B|| and B^ denote the components of magnetic field 
parallel and perpendicular to the mean field Bq, respec- 
tively (see Appendix A of this paper and Appendix C of 
Schekochihin et al. 2009). Here we drop the dissipation 
term for simplicity. 



Finally, from the continuity equation and the assump- 
tion of pressure balance, we can rewrite the above equa- 
tions as 

-Vi X [B.VB||], (18) 



dt 

(see Eqs. 
d b 



dt 



2 -I- PzO_+ Z/t) A-Ken 
and |A14j ). which become 
B 



X [B • VBi] (19) 



dt y/Airnimi 
where 



\/aW± X ( 



^/A-Knimi 



y/ATTTlimi 

(20) 



B = Bo + b, 
b = + 6||Z, 



b = hi 

_ m + Z/r) 



(21) 
(22) 

(23) 

2-f ft(l-HZ/r)- ^^^^ 
Here, Bq (= Bqz) is the mean magnetic field, z is a unit 
vector along the direction of Bq, Z = qi/e (e — \qe\) is 
the ion-to-electron charge ratio, r — Ti/Tf, is the tempera- 
ture ratio, Tie is the average electron number density. The 
vector (in Fourier space) is parallel to z x k^, where 
kj^/Zc^. Therefore, in Fourier space bk lies in the 
plane spanned by z x kj^ and z, which means bk is per- 
pendicular to k^. 

3. EXPECTED SCALING RELATIONS 

In the presence of a strong mean magnetic field, tur- 
bulence energy tends to cascade in the direction perpen- 
dicular to the mean field. As a result, in Fourier space, 
modes with k^ 3> fcy are predominantly excited. In real 
space, characteristic scales parallel to the mean field (Z||) 
tend to be larger than those perpendicular to it (^_l). This 
is referred to as anisotropy of turbulence. 

In the case of standard MHD turbulence, this global 
anisotropy has been studied since early 1980s (Shebalin, 
Matthaeus, & Montgomery 1983). On the other hand, 
Goldreich & Sridhar (1995) showed that there exists a 
regime of turbulence in which a critical balance is main- 
tained between wave motions (with timescale of ~ 
1\\/Va) and hydrodynamic motions (with timescale of 
l±/vi). This is the so-called strong turbulence regime 
and they found a certain relation between k\\ and fcj^, 

2/3 • 

fe|| oc kj^ , in the regime. This scale-dependent anisotropy 
was numerically confirmed by Cho & Vishniac (2000) and 
Maron & Goldreich (2001). Cho & Vishniac (2000) showed 
that this scale-dependent anisotropy can be measured only 
in a local coordinate frame which is aligned with the locally 
averaged magnetic field direction. The necessity of using a 
local frame is due to the fact that eddies are aligned along 
the local mean magnetic field, rather than the global mean 
field Bq. We call this kind of anisotropy as local anisotropy. 

In the case of EMHD turbulence, anisotropy has been 
studied only recently. Dastgeer et al. (2000) and Dastgeer 
& Zank (2003) numerically studied global anisotropy of 2D 
EMHD turbulence. On the other hand, Ng et al. (2003) 
numerically studied local anisotropy of 2D EMHD tur- 
bulence, but used second-order structure functions, the 
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all simulations, so that the dissipation term in the above 
equation is replaced with 7y3(V^)^B. We perform both 
driven and decaying turbulence simulations. 

4.2. The ERMHD Code 

The ERMHD equation is only slightly different from the 
EMHD equation. Therefore, we also adopt a pseudospec- 
tral code to solve the normalized ERMHD equation in a 
rectangular periodic box of size 27r x 27r x Gtt: 



(B • VB) + ry'V^B, 



where 



B = Bo 



b = b| 



(26) 

(27) 
(28) 



limitation of which will be discussed in Section H6.3I In 
CL04, we studied anisotropy of 3D EMHD turbulence. In 
CL04, we for the first time found that critical balance be- 
tween wave motions and hydrodynamic motions also holds 
true in EMHD turbulence and that anisotropy of EMHD, 

1/3 . 

A;|| oc , is stronger than that of standard MHD turbu- 
lence. 

Anisotropy of ERMHD is expected to be similar to that 
of EMHD (Schekochihin et al. 2009). This is under- 
standable from Eq. (I^Dl) . Note that b = b^ + 6||Z and 

b = b^ + ^1/q:6||Z in the equation. The factor a is 
less than or equal to 1. If we assume 6j_ ~ yJTJah\^\, 
then the parallel (or z) component of true magnetic field 
5|| is equal to ^ \/oib± < h±. If a ^ 1 and, hence, 

b ~ b, then Eq. (^0)1 becomes the usual EMHD equa- 
tion. Therefore it is trivial to show that anisotropy of 
ERMHD is similar to that of EMHD. On the other hand, 
if a <C 1 and, hence, b = bj^ -f \/at>|| ^ bx, then the ER- 
MHD equation becomes different from the EMHD equa- 
tion. However, in general, the parallel component of b 
does not contribute much to the term B • Vb, because 
5||Z • V ~ 6||fc|| ^ bj^ • V ^ fo±fc_L in the presence of 
anisotropy, ^ fcj^. This means that what matters for 
energy cascade is the perpendicular component of b, which 
is not directly affected by the value of a. Therefore, even 
in case of a <C 1, we expect that anisotropy of ERMHD is 
similar to that of EMHD. However, this conjecture needs 
to be tested by numerical calculations. 

4. NUMERICAL METHODS 

As we described earlier, EMHD and ERMHD equations 
involve with time evolution of magnetic field only. Since 
magnetic field is divergence-free, we can use incompress- 
ible numerical schemes to solve the equations. 

4.1. The EMHD code 

We adopt a pseudospectral code to solve the normalized 
EMHD equation in a periodic box of size 27r: 

dB 

'dt 

where magnetic field, time, and length are normalized by 
a mean field Bq, the whistler time = L'^{uipe/c)'^ /^e 
(r2e= electron gyro frequency), and a characteristic length 
scale L (see, for example, Galtier & Bhattacharjee 2003). 
The resistivity rj' in equation (j25p is dimensionless. The 
dispersion relation of a whistler waves in this normalized 
units is w = fcfc||i?o- The magnetic field consists of the uni- 
form background field and a ffuctuating field: B = Bo + b. 
The strength of the uniform background field, Bq, is set 
to 1. We use up to 512'^ collocation points. At t = 0, the 
random magnetic field is restricted to the range 2 < fc < 5 
in wavevector space. The amplitudes of the random mag- 
netic field at i = is ~ 1.2. In order to have a more 
extended inertial range, we use hyperdiffusivity for the dif- 
fusion term^. The power of hyperdiffusivity is set to 3 for 

^ We do not observe a strong bottleneck effect, wfiicfi is a common feature in numerical hydrodynamic simulations with hyperviscos- 
ity/hyperdiffusivity. Therefore expect that hyperdiffusivity does not substantially alter physics in the inertial range. However, this should 
be clarified in the future. Discussions on the bottleneck effect for 2-dimensional EMHD can be found in Biskamp, Schwarz, & Celani (1998). 
^ Or, equivalently, we may use fcy = 1,2,3, ... and kj_ = 3,6,9, .... In this case, the size of the computational box is (27r/3) X (27r/3) X {2tt), 
where 2it is the side along the mean magnetic field. 



-V X [(V X B) X B] +T]'\/^B, 



(25) 



i EE y^t, (29) 

where 77' is dimensionless resistivity, b^ and 6||Z are per- 
pendicular and parallel components of magnetic field, re- 
spectively, and definition of a is given in Eq. ([M]) . Note 
that the ERMHD equation is very similar to the EMHD 
equation. Numerical setup is therefore similar to that of 
the EMHD case, except the shape of the simulation box. 

We use an elongated numerical box for ERMHD to sat- 
isfy the condition k± ^ fc||. The numerical box is 3 times 
longer in the direction of the mean magnetic field. The 
wavenumbers along the mean field direction have frac- 
tional values 

= 1/3,2/3,1,4/3,..., (30) 

while those of perpendicular direction have integer value^. 
We perform only decaying turbulence simulations. At t=0, 
Fourier modes with 1/3 < fc|| < 1 and 4.5/-\/2 < fc^ < 

are excited. The strength of the mean magnetic 
field (Bq) is 1 and the rms value of the random magnetic 

field (actually b) at t=0 is ~ 0.071. Therefore, the critical 
balance, bk±/Bok^^ — 1, is roughly satisfied at t=0. For 
comparison, we perform a similar EMHD simulation (Run 
E256D-EL). 

5. RESULTS: EMHD AND ERMHD 

We compare EMHD and ERMHD turbulence in Fig. [TJ 
Left panel of Fig. [T] shows time evolution and energy 
spectra of decaying EMHD and ERMHD runs (see Runs 
E256D-EL, ER256D1-EL, and ER256D8-EL in Table 1). 
The inset shows time evolution of magnetic energy density. 
The solid curve is for ERMHD with ^/a = 1, dotted curve 
for ERMHD with y/a = 1/8, and the dashed curve for 
EMHD. The horizontal axis represents time multiplied by 
^/a and the vertical axis twice the magnetic energy density. 
In case of ERMHD with = 1/8 (i.e. Run ER256D8- 
EL), the vertical axis is not 5^, but P — + b^a. All 
three curves follow a similar decay law. The main plot of 
Fig. [T] shows energy spectra at ^/at = 0.6 (see the arrow 
in the inset). All three spectra are consistent with the 
expected k~'^/^ spectrum. 
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Fig. 1.— Comparison of EMHD (Run E256D-EL) and ERMHD (Runs ER256D1-EL and ER256D8-EL). Left: energy spectra and time 
evolution of magnetic energy density (inset). The energy spectra are taken at the time marked by the arrow in the inset. Note that the 
horizontal axis of the inset is y/at, where « is a parameter that appears in the ERMHD formalism (Eq. 1241 ). The vertical axis of the inset is 
(or, in case of ER256D8-EL, b^). Right: comparison of anisotropy. The calculations are done at the time marked by the arrow in the inset 
of left panel. 



We compare anisotropy of EMHD and ERMHD in right 
panel of Fig. [T] We obtain the plot at ^/ai = 0.6. All three 

curves show similar scaling relations, fc|| cx k^j^"^ . Note 
however that, although we do not show it in this paper, 
anisotropy tends to get stronger as turbulence decays fur- 
ther. We will describe the method of obtaining anisotropy 
and its limitations in Section i j6.1l 

To summarize, we observe that both EMHD and ER- 
MHD turbulence exhibit similar spectra and anisotropy. 
Spectra and time evolution of magnetic energy density of 
ERMHD turbulence seem to be invariant under the fol- 
lowing transformation, 

B ^ B, Vat- (31) 

6. RESULTS: SCALING OF EMHD 

In the previous section, we have studied turbulence in a 
rectangular box elongated along the mean field direction. 
Note that the numerical resolution in the perpendicular 
direction is 256 x 256, which limits the inertial range of 
turbulence. To increase the inertial range we need to in- 
crease numerical resolution in the perpendicular direction. 
Studying ERMHD is more difficult than studying EMHD 
because we need to use an elongated numerical box for 
ERMHD to satisfy the condition k± > k\\. For EMHD, 
there is no necessity of using such an elongated numerical 
box. Therefore, EMHD is more suitable for high resolution 
simulations. Since EMHD turbulence and ERMHD turbu- 
lence exhibit similar scaling relations, we only consider a 
high resolution EMHD simulation in this section. 

We perform a decaying EMHD simulation with 512'^ grid 
points (Run E512D in Table 1). At t=0, Fourier modes 

with 2 < fc < 5 are excited. Here k = \Jk\ + fcy and, 

thus, the initial perturbation is isotropic (i.e. ~ fc^). 
The strength of the mean magnetic field {Bq) is 1 and 
the rms value of the random magnetic field (6) is ~ 1.21. 
Therefore, the critical balance, bk±/Bok\\ = 1, is roughly 
satisfied at t=0. 

In this section, we focus on anisotropy of EMHD turbu- 



lence. CL04 showed that anisotropy of EMHD turbulence 

1 /3 

is consistent with fc|| oc fcy . However, since it is difficult 
to measure anisotropy of EMHD turbulence, more careful 
assessments of anisotropy are necessary. In this subsec- 
tion, we develop and test new techniques for anisotropy. 
We analyze anisotropy of turbulence after turbulence has 
developed a full inertial range. 

The energy spectra in left panel of Figure [2] show how 
the inertial range develops. At t = only large-scale (i.e. 
small fc) Fourier modes are excited. The dashed curve in 
Figure [2] shows the initial spectrum. As the turbulence de- 
cays, the initial energy cascades down to small scales and, 
as a result, small scale (i.e. large fc) modes are excited. 
When the energy reaches the dissipation scale at fc > 100, 
the energy spectrum goes down without changing its slope 
(the dotted and the solid curves). The slope at this stage 
is very close to that of the predicted spectrum: 

E{k) (X fc-^/^ (32) 

For the study of anisotropy, we use the data cube at 
t « 0.46. The solid curve in left panel of Figure rep- 
resents the spectrum at that time. 

The right panel of Figure [2 shows spectra of electric 
field. Both the total electric field (solid curve) and its xy- 
componcnt (dashed curve) show spectra compatible with 
k"^/^. (Note that we assume Bq = Bqz.) This result 
is consistent with the ones obtained by gyrokinetic simu- 
lations (Howes et al. 2008a) and Hall MHD simulations 
(see Dmitruk & Matthaeus 2006; Matthaeus et al. 2008). 
Bale et al. (2005) reported a similar electric fluctuation 
spectrum in the solar wind. 

6.1. Anisotropy: method 1 

In CLG4, we noted that the quantity B^ • Vb; is propor- 
tional to _Bifc||fei, where B^ is the local mean field, b; the 
fluctuating field at scale I, and fc|| the wave number parallel 
to the local mean magnetic field. We obtained B/, by elim- 
inating Fourier modes whose perpendicular wavenumber is 
greater than fc/2. We obtained the fluctuating field b; by 
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Fig. 2. — Spectra of B (left) and E (right) of EMHD turbulence (Run E512D). Left: spectrum at different time points. At t 0, only 
small-wavenumber Fourier modes are excited. At later times, energy gradually cascades down to small scales (i.e. large k regions) and a 
fc— ■^/S inertial range develops. Right: spectrum of electric field E at t=0.46. The solid curve is the spectrum of total electric field and the 
dashed curve is that of electric field perpendicular to the mean field Bp. 




eliminating Fourier modes whose perpendicular wavenum- 
ber is less than fc/2. In CL04, we calculated anisotropy in 
Fourier space: 

1/2 



fe<|k'|<A;+l 



|BL-Vb,|^, 



/c<|k'|<fe+l l^lk' 



(33) 



The reason why we did the calculation in Fourier space is 
that there might be unknown contamination from scales 
other than We plot anisotropy obtained this way in left 
panel of Fig. [31 

We note that, when we obtain the fluctuating field b/ by 
filtering out Fourier modes whose perpendicular wavenum- 
ber is less than k/2 or greater than 2k, we may calculate 
anisotropy directly in real space: 

B..Vb.«W.^fc„«f^Hi_^^V''. (34) 



We plot anisotropy obtained this way in right panel of 
Fig. [31 The result is compatible with the expected 

1/3 

anisotropy, fcii cx kj^ . The limitation of this method is 



that the choice of the filtering wavenumbers, k/2 and 2k, 
is an arbitrary one. When we use filtering wavenumbers 
of fc/\/2 and V2k, we get a slightly weaker anisotropy. 

6.2. Anisotropy: method 2 

In the previous subsection, we assumed that • Vb; (x 
B^k^^bi. Is this really correct? If small-scale eddies are 
aligned with local mean field, the above assumption should 
be true. However, if the major axes of small-scale eddies 
are substantially misaligned with local large-scale mean 
field, the quantity Bl • Vb/ may sample perpendicular 
wavenumber of the eddies. In this subsection, we investi- 
gate alignment of small-scale eddies with respect to local 
mean magnetic field. 

We first visualize the alignment effect. Left panel of 
Fig. m shows a snapshot of magnetic field structure in 
a plane parallel to the global mean magnetic field. The 
global mean magnetic field is parallel to the horizontal 
axis in the Figure. The white lines in Figure represent 
local mean magnetic field. We obtain the local mean field 
by filtering out Fourier modes with k > 15. The contours 
represent magnetic energy density of small-scale field. We 
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Fig. 4. — Anisotropy based on local correlation lengths (Run E512D). Left: visualization. White lines denote local mean magnetic field 
obtained by filtering out large-fc Fourier modes. Contours show structure of small-scale eddies obtained by retaining Fourier modes near 
the scale of interest and filtering out all other scale Fourier modes. The global mean field Bq is parallel to the horizontal axis. Right: a new 
method based on correlation lengths in local frame. The parallel correlation length iy, for example, is the average of the correlation length 
measured along the local mean magnetic field direction at each point. (Note that the local mean magnetic field directions are marked by 
white lines in the left panel). In Run E512D, A^side = 512. 



obtain the small-scale field by filtering out Fourier modes 
with /c < 15 or fc > 60. Since the alignment effect, if 
any, should be 3-dimensional effect in nature, we may not 
be able to visualize the effect correctly in a 2-dimensional 
plane. Nevertheless, we can observe that small-scale ed- 
dies are well aligned with the local mean field directions. 

If small-scale eddies are well aligned with the local mean 
field, the correlation length along the local mean field 
should be larger than that of perpendicular directions. 
Furthermore, when the alignment effect is present, the 

parallel correlation length should follow a l^"^ (or, k^^^"^) 
power law. We investigate the behavior of parallel and 
perpendicular correlation lengths by changing the scale /0 
We define the correlation lengths as the distances at which 
the 2-point correlation function (b;(ri) • b/(r2)) drops by 
50%. The parallel correlation length is measured along 
the directions of local mean field, B^, and the perpen- 
dicular one for perpendicular directions. In right panel 
of Fig. 21 we plot the parallel and perpendicular correla- 
tion lengths as functions of small-scale eddy size I cx 1/k. 
The parallel correlation lengths (upper curve) is of course 
longer than the perpendicular ones (lower curve). On av- 
erage, the parallel correlation lengths seem to follow the 
expected l^^"^ scaling. However, the slope is shallower than 
1 /3 for small-size eddies and steeper than 1/3 for large-size 
eddies. Roughly speaking, the perpendicular correlation 
lengths follow scaling, which is reasonable. 

All in all, although it may not be not a perfect method, 
the filtering method reasonably reveals anisotropy of 
EMHD turbulence. The anisotropy of EMHD turbulence 
revealed by the filtering process is consistent with the 



6.3. Anisotropy: Multi-point second-order structure 
functions 

We may visualize scale-dependent anisotropy using the 
second-order structure function calculated in the local 
frame, which is aligned with local mean magnetic field Bl: 

SF2(r||,r_L) =< |B(x + r)-B(x)|2 > 

avg. over x; 

(35) 

where r = 7'||f || + r±rj_. The vectors fy and f are unit 
vectors parallel and perpendicular to the local mean field 
Bl, respectively. See Cho & Vishniac (2000) for the de- 
tailed discussion of the local frame. Left panel of Fig. [5] is 
an example of such visualizatiorQ. The horizontal axis cor- 
responds to the local mean field directions. The shapes of 
contours illustrate scale-dependent anisotropy: the large 
eddies are roughly isotropic and smaller eddies are more 
elongated. However, we should be careful when we de- 
rive a quantitative scaling relation for anisotropy from the 
contour diagram. 

The 2-point second-order structure function of a vari- 
able A 



SF2{r) = {8Arf =< |^(x + r) - A(x)|2 > 



(36) 



is in general related to the energy spectrum of the variable, 
EA{k), through 

SF2{r) = (SAr)^ ^ kEAik), (37) 
where k cx 1/r is the wavenumber. Therefore, when 
EA{k) cx fc~™, we have 

S'F2(r) cxr'"~^ (38) 
For example, in fully developed Kolmogorov turbulence 
{E{k) cx fc~^/'^), the scaling relation of the second-order 
longitudinal structure function is given by 

SF2{r) = {5vrY - kk""^'^ cx fc^^/s ^ ^2/3^ (39) 
However, when the slope of the turbulence spectrum is 
steeper than k~^ , the relation in Eq. ([57|) becomes invalid 

^ We use a similar filtering method as in the previous subsection. Wc obtain the local mean field, B^, by eliminating Fourier modes whose per- 
pendicular wavenumber is greater than k/2. We obtain the fluctuating field, bj, by eliminating Fourier modes whose perpendicular wavenumber 
is less than k/2 or greater than 2k. Note that I = 27r/A;, or in terms of grid units I = N^i^^/k, where N^i^^ = 512 in Run E512D. 

^ Note however that we do not use the usual 2-point second-order structure function for the plot. Instead, we use the 5-point second-order 
structure function that will be discussed later in this subsection. 



CX kY^ scaling. 
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Fig. 5. — Anisotropy from the second-order structure function (Run E512D). Left: contour diagram obtained from the 5-point second-order 
structure function. Middle: comparison of the 3-point SF2 and the 5-point SF2. The parallel structure function SF2{r^^,0), for example, can 
be obtained from values of SF2 on the horizontal axis of the contour diagram. Right: the relation between the semi-minor axis (oc l/fcj^) and 
the semi-major axis (oc l/fcy) of the contours. 



and SF2 (r) cx regardless of the slope of the energy spec- 
trum (see Appendix B). In EMHD, the energy spectrum 

expressed in terms of k± scales as E{k±) cx kj_ and we 
expect that the second-order structure function scales as 
SF2iO,r^)<xr^/\ (40) 
On the other hand, the energy spectrum expressed in 
terms of /cj| is expected to be E{k^^ ) oc k'^^ when anisotropy 

1 /3 

scales as fcy oc A;/ (see CL04). Therefore, the one-to-one 
correspondence between the second-order structure func- 
tion and the energy spectrum is not valid for parallel direc- 
tion. This means we cannot use the second-order structure 
function to reveal true anisotropy. 

We have shown that the 2-point second-order structure 
function is not suitable for quantitative study of anisotropy 
in EMHD turbulence. However, it is possible to construct 
multi-point second-order structure functions that can be 
used for variables with steep energy spectra (see Falcon, 
Fauve, & Laroche 2007; Lazarian & Pogosyan 2008; see 
also Appendix C). The 3-point SF2 was used by Falcon et 
al. (2007) and Lazarian & Pogosyan (2008). The 3-point 
SF2 will work for energy spectrum as steep as ~ k~^. 
In Appendix C, we discuss how to construct 4-point and 
5-point structure functions. The 5-point SF2 works for 
energy spectrum as steep as ~ k~^. 

Left panel of Fig. O is obtained with this 5-point struc- 
ture function. Middle panel of the Figure shows the 
parallel structure functions, SF2{r^^,0), and the perpen- 
dicular structure functions, SF2{0,r±). We plot the re- 
sults of 3-point (dotted lines) and 5-point (solid lines) 
structure functions. In the perpendicular direction, both 
structure functions reasonably follow the expected scaling 
of . In the parallel directions, however, we observe 
that, when r|| > 10, the structure functions are much 
shallower than the expected scaling of r|, which is from 
5^2 ^11 -£'(^11) oc oc r|. In the parallel directions, the 
5-point structure function is steeper than the 3-point one. 

One should note that it is very difficult to obtain a well- 
defined power-law scaling for the parallel direction. This is 
because the inertial range in the parallel direction, if any, 
is extremely short. In the perpendicular direction, the in- 
ertial range spans from the outer scale {k± ~ 3) to the 

1/3 

dissipation scale {k± ~ 100). Suppose that fcy oc k^ is 
the true anisotropy. Then, in the parallel direction, the in- 
ertial range spans from fc|| ~ 3 to ~ 3 x (100/3)1/3 ^ 10. 



Therefore, it is virtually hopeless to reveal true anisotropy 
from contour diagram. Indeed, right panel of Fig. [5] shows 
that the anisotropy derived from the relation between 
semi-major axis and semi-minor axis of contours is not 

1 /3 

really consistent with the fc|| oc kj^ relation. Both the 
solid curve (5-point) and the dotted curve (3-point) show a 
similar scaling. The average slope is approximately ~ 0.5. 
Therefore, we can conclude that true anisotropy is similar 
to or stronger than this. Note that the dashed curve (2- 
point) show a milder anisotropy (or, a steeper slope). This 
means that multi-point structure functions are indeed bet- 
ter for steep spectrum. Although structure functions may 
not be suitable to reveal true anisotropy for EMHD, it 
is promising that multi-point structure functions seem to 
resolve steeper spectrum better. 

6.4. Decaying vs. Forced EMHD turbulence 

It is generally true that both decaying and driven tur- 
bulence show a similar scaling. However, in the presence 
of strong mean magnetic field, it is not clear whether or 
not they really show a similar scaling. Consider a de- 
caying strongly magnetized EMHD turbulence. Let us 
assume that critical balance is roughly satisfied at t=0. 
The strength of random magnetic field drops as time goes 
on. But, the strength of the mean magnetic field remains 
same. Therefore, critical balance will be soon destroyed. 
Of course, there is a narrow window of time during which 
critical balance is satisfied and we can study critically bal- 
anced EMHD turbulence during the time interval. Never- 
theless, it is worth comparing driven turbulence and de- 
caying turbulence. In this section, we show that forced 
EMHD turbulence exhibits similar scaling relations as de- 
caying EMHD turbulence. 

Numerical setup for driven EMHD turbulence is similar 
to that of decaying EMHD turbulence. Forcing is done 
in Fourier space. The forcing term consists of 21 Fourier 
components with 2 < k < vT2. The peak of energy in- 
jection is at k « 2.5. The forcing is statistically isotropic. 
The strength of the mean magnetic field Bq is 1. We ad- 
justed the amplitudes of the forcing components so that 
6 w 1. Therefore, the driven EMHD turbulence satisfies 
the condition for critical balance. This run is designated 
as Run E256F in Table 1. 

Inset in Fig. [6] shows time evolution of magnetic energy 
density (in fact = Bq + b^). The magnetic energy 
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Fig. 6. — Driven EMHD turbulence (Run E256F). The strength of the (global) mean field is 1 (i.e. Bq = 1). The strength of the fluctuating 
magnetic field (b) is maintained to be ~ 1. Since driving is isotropic, the critical balance is satisfied. Left: time evolution of magnetic energy 
density (inset) and energy spectra. The spectra are taken at 3 different time points marked by the arrows in the inset. Right: anisotropy. 
Line styles of the curves are the same as those of the arrows in the inset. 



reaches a stationary state after t ^ 1. The main plot of 
Fig. [S] shows magnetic energy spectrum at 3 different time 
points marked by the arrows in the inset. Line styles of the 
spectra are the same as those of the arrows. The slopes of 
all 3 energy spectra are consistent with —7/3. 

We plot anisotropy in right panel of Fig. [6l We use the 
technique described in Section ^ 16. II Anisotropy of driven 

1 /3 

EMHD turbulence is also consistent with the fcy oc kj^ 
scaling. All in all, driven MHD turbulence shows similar 
spectrum and anisotropy as decaying EMHD turbulence. 

6.5. Inverse cascade of EMHD turbulence 

In Wareing & HoUerbach (2009) the inverse cascade was 
reported as a part of decaying 2D EMHD turbulence. The 
energy spectrum in Shaikh & Zank (2005) also shows a 
clear evidence of inverse cascade in driven 2D EMHD tur- 
bulence. However, the difference of 3D and 2D turbulence 
makes one wonder whether the inverse cascade is present 
in decaying 3D EMHD turbulence. 

We perform another decaying EMHD turbulence simu- 
lation. The numerical setup is almost identical to the one 
described in Section §4.1. However, we use a different ini- 
tial condition. At t = 0, Fourier modes with 16 < A: < 32 
are excited (see the solid curve in Left panel of Fig. [7]). 
The numerical resolution is 256"^. This run is not listed in 
Table 1. 

Left panel of Fig. [7] shows spectral behavior of the de- 
caying EMHD turbulence. The dotted line shows the spec- 
trum shortly after the simulation. The dotted line clearly 
shows that most of the energy cascades down to small 
scales. Of course, the energy cascading down to small 
scales will dissipate away below fc > 80. It also shows 
that some of the energy goes to larger scales, the amount 
of which is smaller than that in the 2D EMHD case (see 
Wareing & HoUerbach 2009). At later times, the peak to 
the energy spectrum moves to larger scales, which prob- 
ably means that we see a self-similar decaying of energy 
that leads to the increase of the integral scale (see, for 
example, Biskamp 2003), rather than a signature of the 
inverse energy cascade that is observed in 2D hydrody- 



namic turbulence. Nevertheless, we clearly observe that 
EMHD turbulence can generate more coherent magnetic 
field from a smaller scale, hence less coherent, magnetic 
field. Although this phenomenon is not the inverse cas- 
cade per se, we can still call it inverse cascade because 
it does show that small amount of energy goes to larger 
scales. However, in order to avid confusions, we will use 
the term 'small amount of inverse cascade' whenever pos- 
sible. 

We compare our results with the spectral behavior of 
the decaying 3D MHD turbulence in Right panel of Fig. [71 
At t = velocity field has a flat spectrum between k = 16 
and A: = 32 (solid curve). The magnetic field at t = 
has only the uniform component, the strength of which is 
set to 1. The numerical resolution is also 256'^ and this 
run is not listed in Table 1. At later times, the velocity 
field generates fluctuating magnetic field mainly between 
fc = 16 and k = 32. We see that MHD turbulence also 
shows small amount of inverse cascade. Overall spectral 
behavior of decaying MHD turbulence is similar to that 
of decaying EMHD turbulence. However, there are also 
differences. For example, the gradual shift of the peak of 
the energy spectrum to larger scale is less pronounced in 
the MHD case. The slopes of the energy spectra on large 
scales are also different. 

7. HIGH-ORDER STATISTICS AND BISPECTRUM OF EMHD 

In this section, we present other statistical properties of 
EMHD turbulence. We use a data cube from Run E512D, 
which is the the same data cube as we considered in Sec- 
tion ^ 

7.1. High-order structure functions 

High-order structure functions are used for the study of 
intermittency, which refers to the non-uniform distribu- 
tion of structures. The structure functions of order p for 
magnetic field is defined by 

SFp{r) =< |B(x)-B(x + r)|P > 

avg.over x • 

(41) 

Traditionally, researchers use high-order structure func- 
tions of velocity to probe dissipation structures of tur- 
bulence. In fully developed hydrodynamic turbulence, 
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Fig. 7. — Spectra of decaying turbulence. The mean field {Bq) is set to 1 in both cases. Left: EMHD turbulence. At t = 0, the magnetic 
spectrum is flat for 16 < A; < 32 (solid curve). At later times, energy cascades down to smaller scales. At the same time, some energy cascades 
inversely up to larger scales. Right: MHD turbulence. At t = 0, the velocity spectrum is flat for 16 < fe < 32 (solid curve). The magnetic 
field has only the uniform component (So) at t = 0. At later times, random magnetic field is generated and most of the turbulence energy 
cascades down to smaller scales. We also observe inverse cascade. Note that, in the right panel, the solid line is for velocity at t=0 and other 
lines are for magnetic field at later times. 




Fig. 8. — Intermittency (Run E512D). Scaling exponents from the 3-point, the 4-point, and the 5-point structure functions follow the 
She-Leveque scaling. However, those from the standard 2-point structure function show a different scaling. The <^p's shown here are the 
relative scaling exponents, Cp/Cs- 



the (longitudinal) velocity structure functions SFp 
([v(x + r) — v(x)] • f)P >=< (Svr)'^ > are expected to 
scale as r^p . One of the key issues in this field is the func- 
tional form of the scaling exponents ^p. There are several 
models for Cp- Roughly speaking, the dimensionality of 
the dissipation structures plays an important role. 

Assuming 1-dimensional worm-like dissipation struc- 
tures, She & Leveque (1994) proposed the scaling relation 

Cp^^=p/9 + 2[l-(2/3)^'/3] (42) 
for incompressible hydrodynamic turbulence. On the 
other hand, assuming 2-dimensional sheet-like dissipation 
structures, Miiller & Biskamp (2000) proposed the relation 

Cp*'^ =p/9 + l- (l/3)f/3 (43) 

* In hydrodynamic turbulence, the relative scaling exponents are expressed relative to since Kolmogorov's four-fifth law is valid for the 
third-order longitudinal structure function: SF3 oc r, where the constant of proportionality is (-4/5) times the energy injection rate. Therefore, 
it is natural to use (^3 for the relative scaling exponents. In MHD, there also exists a similar exact relation for a third-order structure function 
expressed in terms of Elsasser variables (Politano & Pouquet 1998). However, according to the exact correlation law obtained by Galtier 
(2008), there is no simple relation between SF3 and r in EMHD and, therefore, the use of 1^3 for the relative scaling exponents is not justified. 
Nevertheless, for the sake of comparison with existing models, we use (^3 for the Figure. 



for incompressible magneto-hydrodynamic turbulence. 
There are other models for intermittency. But, in this 
paper, we only consider above mentioned two models be- 
cause they are relevant to incompressible variables. 

In this paper, we only consider structure functions in 
perpendicular directions. That is, in our calculations, the 
vector r (see Eq. [IT]) is perpendicular to the direction of 
the local mean field, B l ■ The reason why we do not con- 
sider parallel direction is that structure functions are not 
suitable for revealing true scaling relation in parallel direc- 
tions. In Fig. [5] we plot the relative scaling exponent Cp/Cs 
for EMHD turbulenc^. It is interesting that the 2-point 
structure functions show a different scaling exponents com- 
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pared with other multi-point structure functions. This re- 
sult is not surprising because we already observed that the 
2-point second-order structure function shows a different 
behavior compared with other multi-point structure func- 
tions in right panel of Fig. O The scaling exponents based 
on the 3-point, the 4-point, and the 5-point structure func- 
tions are consistent with the She-Leveque model. 

7.2. Bispectrum 

In astronomy, the bispectrum is a tool widely used in 
cosmology and gravitational wave studies (Fry 1998; Scoc- 
cimarro 2000; Liguori et al. 2006). Here, we briefly de- 
scribe the definition of the bispectrum (see Burkhart et 
al. 2009). The bispectrum is closely related to the power 
spectrum. The Fourier transform of the second-order cu- 
mulant, i.e. the autocorrelation function, is the power 
spectrum while the Fourier transform of the third order 
cumulant is known as the bispectrum. In a discrete sys- 
tem, the bispectrum is defined as: 

BS{kuk2)^ Mki)-A{k2)-A*{ki + k2) 



ki = 



1st ko 



(44) 

where ki and ^2 are the wave numbers of two interacting 
waves, and A{k) is the original discrete data with finite 
number of elements with A*{k) representing the complex 
conjugate of A{k). As is shown in Equation (UH), the bis- 
pectrum is a complex quantity which will measure both 
phase and magnitude information between different wave 
modes. 

Original formalism for bispectrum is suitable for scalar 
variables. Since we deal with magnetic field, which is a 
vector quantity, we simply take z-component of magnetic 
field. Note that Bq — BqZ. In Fig. [51 we plot bispec- 
trum of driven EMHD (left panel), decaying EMHD (mid- 
dle panel), and driven MHD (right panel). For the driven 
ERMHD and MHD turbulence cases, we take the data 
after turbulence has reached statistically stationary state 
(i ~ 3 and t ^ ib, respectively). The decaying EMHD 
case was taken at i ~ 3 (see the inset of Fig. [1]). Both 
driven EMHD (left panel) and decaying EMHD (middle 
panel) cases show similar bispectrum. However, driven 
MHD (right panel) shows a substantially different bispec- 
trum. All cases show that amplitudes of bispectra are 
highest when ki — k2, which means a high correlation of 
modes at fci — k2- The amplitudes of bispectra, hence 
strengths of nonlinear interactions, drop as we move away 
from the diagonal lines (i.e. the lines of ki — k2). The 
shape of isocontours depends on the nature of nonlinear 
interactions. The case of the standard MHD (right panel) 
shows wider isocontours than EMHD cases, which means 
nonlinear interactions in standard MHD turbulence drop 
more slowly as we move away from the ki = k2 line. This 
demonstrates that nonlinear interactions in ERMHD and 
MHD cases are very different. 

8. COMPARISON WITH THE SOLAR WIND TURBULENCE 

In this subsection, we compare statistics of 
MHD/EMHD turbulence and turbulence in the solar wind. 
We consider quantities related to the probability distribu- 
tion functions (PDFs) of fluctuations in increments of the 



magnetic field strength B 

dB= (B(x + r) -B(x)),/(S(x))„ (45) 
where B = |B| is the strength of magnetic field and the 
angled brackets (...)a: denote average taken over x. In this 
section, we are only concerned with characteristics of tur- 
bulence directly measured in space plasmas. 

8.1. PDFs for dB 

In the solar wind, the PDFs for increments of the mag- 
netic field strength, dB, is measured by 

dBir) = {Bit + t)- B{t))t/{B{t))u (46) 
where B{t) is the strength (or the averaged strength) of 
magnetic field at time t and (...)( denotes average taken 
over time. The observed PDFs for dB on scales from 1 
hour to 128 days are well described by the Tsallis distribu- 
tion (Burlaga, F-Vinas, & Wang 2007; Burlaga & F-Vinas 
2005), which is given by 



y{x) ^ A 



'1/(9-1) 



(47) 



where A, q, and w are constants. The function is propor- 
tional to a Gaussian function for small x, 



2 \ 2 



as X 



0, 



(48) 



(49) 



y{x)K Ail - — ] wAexp( 

and a power law for large x, 

y[x) lyi x~'^/'^'^~^\ as X — 
The shape of the function is determined by the values of 
q and w. In the limit of g ^ 1, the function reduces to a 
Gaussian distribution. Therefore the value of g is a mea- 
sure of non-Gaussianity. The value of w is related to the 
width of the distribution. Since r > 1 hour, the measured 
PDFs in Burlaga et al. (2007) and Burlaga & F-Vinas 
(2005) refiect fiuctuations of MHD turbulence. 

Fig. [To] shows PDFs of MHD (top panel) and EMHD 
(middle and bottom panels) turbulence. The solid lines 
are PDFs of dB in our numerical data and the dotted 
curves are fits of the PDFs to the Tsallis distribution. We 
use the Levenberg-Marquardt algorithm (Levenberg 1944; 
Marquardt 1963; see Press et al. 1992) for fitting. Each 
curve in the panels represents the PDF for a particular 
separation. As we move from the lowest curve in each 
panel, the separations in grid units are 2, 6, 10, 17, 29, 50, 
88, respectively. Each curve is displaced vertically from 
the one below by a factor of 10. We set the value of each 
PDF at dB = to 1 before displacement, for the sake of 
clarity. 

The PDFs may depend on the direction of r in Eq. (145]) . 
That is, the PDFs for the parallel direction can be differ- 
ent from those for perpendicular direction. Therefore, we 
calculate PDFs for parallel and perpendicular directions 
separately. We try both global and local frames to define 
parallel and perpendicular directions. Overall, we consider 
4 cases: parallel direction in global frame (r||i3o), parallel 
direction in local frame (r||_BL), perpendicular direction 
in global frame (r _L -Bq); and perpendicular direction in 
local frame (r _L Bl). 

Our calculations show that there is no big difference in 
the PDFs of the above mentioned 4 cases. General trend 
is that the PDF is close to a Gaussian function when sepa- 
ration is large (upper curves in each panel) and it deviates 
from a Gaussian function when the separation is small 
(lower curves in each panel). We note however that, for a 
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Fig. 9. — Bispectrum of EMHD (left and middle panels) and standard MHD (right panel). Driven EMHD turbulence at t ~ 3 (left panel; 
Run E256F) and decaying EMHD turbulence at t ~ 3 (middle panel; Run E256D) have similar bispectra. However, driven MHD turbulence 
at t ~ 45 (right panel; Run MHD256F) shows a different bispectrum. 
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Fig. 10. — The PDFs for increments of magnetic field and Tsallis fit. The Tsallis distribution fits the PDFs well. The PDF is close to 
a Gaussian distribution when separation r is large (upper curves in each panel) and it deviates from a Gaussian function as the separation 
gets smaller (lower curves). Runs E512D (middle and lower panels) and MHD512F (upper panel) are used. Note that Bq denotes the global 
mean field and Bj^ the local mean field. In each panel, curves correspond to r = 2 (lowermost curve), 6, 10, 17, 29, 50, 88 (uppermost curve), 
respectively. 



given separation (for example, see the lowest curve in each 
panel, which corresponds to r = 2 grid units), the widths 
of the PDFs vary, depending on the direction of r. 

We can confirm quantitatively the trend that the PDF 
is close to a Gaussian distribution when separation is large 
and it deviates from a Gaussian distribution when the sep- 
aration is small. In Fig. [11] we plot the values of q. The left 
panel is for EMHD turbulence (Run E512D) and the right 
panel for MHD turbulence (Run MHD512F). As we men- 
tioned, the value of q is close to 1 when the PDF is close 
to a Gaussian distribution. Indeed, when the separation is 
large, q is very close to 1. The value of q deviates from 1 
as the separation gets smaller, which means that the PDF 
deviates from a Gaussian distribution. The overall trend 
is consistent with the observations of the MHD-scale fluc- 
tuations in the solar wind (see Burlaga & F-Vinas 2005; 
Burlaga et al. 2007). But, this should be taken as a very 



approximate statement, because the behavior of q in the 
solar wind is very complicated. In each panel of Fig. 111! we 
plot the q values of 4 different directions mentioned above. 
In general, all 4 cases show a similar trend, although the q 
values for the parallel cases are slightly larger than those 
for perpendicular cases. 

In Fig. [T^l we plot the values of w. The left panel is for 
EMHD turbulence (Run E512D) and the right panel for 
MHD turbulence (Run MHD512F). It is interesting that 
the values of w in EMHD cases show steeper dependence 
on separation than those of MHD cases. One should also 
note that the behavior of q as well as w for parallel di- 
rections in EMHD cases may be dominated by large-scale 
fluctuations, because the energy spectrum E{k^^) is too 
steep (see Section ^6.31 and Appendix B). 

Fig- HH shows that the values of w for perpendicular di- 
rection (dotted lines) are systematically larger than those 
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separation (grid units) separation (grid units) 

Fig. 11. — The q parameter of Tsallis distribution. When the separation is large the values of q are very close to 1, while they are larger 
for smaller separations. Note that q = 1 corresponds to a Gaussian distribution. Runs E512D (left panel) and MHD512F (right panel) are 
used. 
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Fig. 13. — Angular dependence of the w parameter of Tsallis distribution. Left: The PDFs for dB and Tsallis fit. The curves in each 
panel correspond to 9 = 0° (lowermost curve), 15° , 30° , 45° , 60° , 75° , and 90° (uppermost curve), respectively. 8 is the angle between local 
mean magnetic field and the direction of r. Right: Dependence of w on 9. The horizontal axis is 8 in degrees. Runs E512D (EMHD) and 
MHD512F (MHD) are used. All calculations are done in local frame. 



for parallel directfon (solid lines). Fig. [T3l shows that w parameter indeed depends on the direction of r. The 
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curves in the left panel are the PDFs for different values 
of 6, which is the angle between the local mean magnetic 
field and the separation vector r. The curves in top or 
bottom panel correspond to 6* = 0° (lowermost curve), 
15°, 30°, 45°, 60°, 75°, and 90° (uppermost curve), respec- 
tively. We can clearly see that the width of the PDF gets 
larger as the angle increases. Note that the w parameter 
of Tsallis distribution is related to the width of the PDF. 
The right panel quantitatively shows dependence of w on 
6. The w parameter is smallest when 6 — Q (^parallel 
direction) and it is largest when 9 = 90° (perpendicular 
direction). 

8.2. PDFs for wavelet transform coefficients 

The wavelet transform is sometimes used for the study 
of the magnetic field fluctuation in the solar wind. For ex- 
ample, Alexandrova et al. (2008) used the Morlet wavelet 
transform 

W,(r,i)= ^S,(i,)^[(t,-i)/r], (50) 
j=o 

where Bi(tj) is the ith component of the magnetic field 
at tj = to + j X At, and ip{u) oc cos(6u) exp(— u^/2) is 
the Morlet wavelet. In this subsection, we briefly consider 
PDFs of the Morlet wavelet transform coefficient Wi- 

We perform the Morlet wavelet transform using our nu- 
merical data. We only consider the parallel (or z) compo- 
nent of magnetic field in global frame. Note that, in our 
numerical calculations, r in Eq. ()50|) is proportional to sep- 
aration r and the global mean magnetic field Bq is parallel 
to z-direction. Fig. [HI shows that the Tsallis distributions 
roughly fit the PDFs of EMHD turbulence (Run E512D). 

The horizontal axis is W^(T)/(Wz(r)2)j/^ In general the 
PDFs for large separations (upper curves in each panel) 
are close to Gaussian distributions and those for small 
separations (lower curves in each panel) show departure 
from Gaussian ones. The PDFs for the parallel direction 
(lower panel) show stronger departure from Gaussian dis- 
tributions, which is confirmed by the large values of q for 
small separations (not shown in this paper). In the right 
panel of the Figure, we show the flatness (or Kurtosis) as 
a function of the scale, which is defined by 

{wt)/{w!r. (51) 

We show both MHD and EMHD cases in the same plot. 
In the case EMHD, we show the flatness of both dB^ and 
Wz In the case of MHD, we show the flatness of dB^ only. 
When the separation is large, the flatness is very close 
to the Gaussian value of 3. In general, it increases as the 
separation decreases, which is consistent with observations 
(Carbone et al. 2004; Alexandrova et al. 2008; see also 
Burlaga & F-Vinas 2005). Note that flatness of of 
EMHD turbulence in the parallel direction shows a strong 
dependence on the separation. 

9. DISCUSSIONS 

9.1. Anisotropy of EMHD 

Anisotropy of EMHD is impossible to measure using the 
traditional second order structure function. In CL04, we 

^ Note however that we considered turbulence in strongly magnetized 
to the local mean field B/,. 



proposed to a different technique which allowed us to ac- 

1 /3 . 

tually get /c|| ~ fcj^ , which conflrmed our theoretical pre- 
diction. 

However, as the problem of measuring anisotropy is im- 
portant, in the present paper we provided not only higher 
resolution simulations, but also a few new techniques of 
measuring anisotropy. All these techniques provided re- 
sults consistent with our earlier finding. 

Unlike CL04, in this paper we did not limit ourselves by 
the simulations of the decaying turbulence, but also stud- 
ied driven turbulence. We did not observe differences in 
the slope or spectrum between the two cases. 

9.2. EMHD and Kinetic Alfven wave turbulence 

As we have shown in §2, there are two approaches to 
describing the turbulence over scales less than the proton 
gyroscale. The traditional approach assumes that the term 
proportional to V • v (see Eq. [T3]) is negligible, while the 
approach in Schekochihin et al. (2009) insists on keeping 
this term. We performed numerical simulations with and 
without the term and obtained virtually identical results 
for the turbulence spectrum and anisotropy. This result 
shows importance of EMHD to analyze coUisionless plas- 
mas. 

CL04 derived the anisotropy of the EMHD turbulence 
and confirmed the obtained scaling numerically. Recently, 
Howes et al. (2008a) used a gyrokinetic code and ob- 
tained the expected fc"''/'^ energy spectrum for magnetic 
fiuctuations below the proton gyroscale (see also discus- 
sions in Matthaeus et al. 2008 and Howes et al. 2008b). 
Anisotropy has not been studied with a gyrokinetic code 
yet. 

9.3. MHD and EMHD turbulence 

Within the paper we have employed a set of different sta- 
tistical measures to characterize EMHD turbulence. We 
compared the results with those obtained for MHD turbu- 
lence. We found both similarities and differences between 
the two types of turbulence. 

MHD turbulence and EMHD turbulence have different 
energy spectra and anisotropy in spite of the fact that 
their scaling relations are derived from the same princi- 
ples: constancy of energy cascade and critical balance. 
Our confirmation of the anisotropy scaling predicted in 
CL04 testifies that the importance of the critical balance 
goes beyond MHD. 

Surprisingly, the high-order statistics shows a similar 
scaling: scaling exponents of both MHD turbulence (see 
Cho, Lazarian & Vishniac 2002) and EMHD turbulence 
are compatible with those of the She-Leveque scaling. 
Does it mean that the She-Leveque scaling is so univer- 
sal that it fails to distinguish different types of turbulence? 
What is the underlying reason for this universality? These 
are still open questions to be addressed by the future re- 
search. 

Interestingly enough, the PDFs for increments of mag- 
netic field are well described by the Tsallis function. 
The correspondence between the space measurements and 
Tsallis description can be found, for example, in Burlaga 
& F.-Vinas (2005), Leubner & Voros (2005), and Burlaga 

iium and that the SFp we measured are for directions perpendicular 
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Fig. 14. — The PDF and flatness (Kurtosis) of the Morlct wavelet transformation coefficient VV^, where z is parallel to the global mean 
field Bq. Left: the Tsallis distribution approximately flts the PDFs. In each panel, curves correspond to r = 2 (lowermost curve), 6, 10, 17, 
29, 50, 88 (uppermost curve), respectively. Right: flatness is larger for smaller separation. When the separations are large, the measured 
values of the flatness are very close to 3, which is the same as the flatness of the normal distribution. 



et al. (2007). Now wc confirmed the correspondence with 
direct numerical simulations. 

Another statistics, namely, the bispectrum is very differ- 
ent for MHD and EMHD. The bispectrum measures non- 
linear interactions. The difference observed confirms the 
difference in cascading of the two types of turbulence. 

9.4. Inverse cascade in EMHD 

Our simulations in §6.5 show the existence of small 
amount of inverse cascade in decaying EMHD turbulence. 
Although we have not performed simulations with driven 
turbulence, we expect that small amount of inverse cas- 
cade is an intrinsic part of the EMHD cascade. That is, 
we expect that generation of small amount of larger scale, 
hence more coherent, magnetic field from small scale tur- 
bulence is an intrinsic part of the EMHD turbulence. 

This kind of inverse cascade has important implications 
in the 3D MHD case. When there exists turbulent veloc- 
ity cascade, this kind of inverse cascade enables generation 
of coherent field on the outer scale of turbulence from a 
small-scale seed magnetic field (see Cho et al. 2009). Note, 
however, that the existence of turbulent velocity cascade 
that can amplify the magnetic field is assumed here. In 
EMHD, on the contrary, the fluctuations of magnetic field 
and velocity are connected with each other from the very 
beginning. Nevertheless, one can speculate that the small 
amount of inverse cascade in EMHD turbulence can help 
to create large scale magnetic structures when the driving 
is at small scales. An interesting implication of this kind 
of inverse EMHD cascade could be its contribution to the 
creation of the dipole field, when magnetic field is gener- 
ated by the thermomagnetic instability in the crust of a 
neutron star. 

9.5. Comparison with observations 

Space plasma measurements allow for a direct compari- 
son with the results of turbulence simulations. This com- 
parison is essential, as numerically turbulence can be stud- 
ied only for relatively small Reynolds and Lundquist num- 
bers. As a result, validation of the obtained scaling in 



realistic astrophysical environments is absolutely (essential 
(see an extended discussion of the issue in Lazaiian et al. 
2009). 

It is encouraging that the spectra and the PDFs of the 
increments of magnetic field strength obtained numerically 

show correspondence with observations. We believe that 
more detailed comparisons are necessary. 

9.6. EMHD turbulence anisotropy and the ADAF model 

The issue of why accretion disks around black holes are 
not as luminous as one would expect is a burning ques- 
tion, addressing of which is required for explaining the 
low luminosity of black hole environments in the centers 
of galaxies. One of the ideas proposed was that of Advec- 
tion Dominated Accretion Flows (ADAFs) (Narayan & Yi 
1995; Narayan et al. 1998). According to this idea the low 
luminosity of the accreting material is due to the low rate 
of the transfer of energy of the turbulent accreting flow to 
electrons. It is postulated in the model that protons carry 
the lion's share of the flow energy and thus the emission 
is suppressed. Is it so? 

Quataert & Gruzinov (1999) discussed transition from 
standard Alfvenic MHD turbulence to EMHD turbulence 
in advection dominated accretion flows. As Alfvenic tur- 
bulence reaches the proton gyroradius scale, some of the 
turbulence energy goes to protons through coUisionless 
damping. The major heating mechanism for protons right 
above the proton gyroscale is the transit time damping 
(TTD) caused by non-zero parallel magnetic field fluctu- 
ations. Heating of protons is sensitive to Pi w (wj/ii^)^, 
where va is the Alfven speed. When /3i » 1, more protons 
are available for efficient interaction with Alfven waves. 
Therefore, most of the turbulence energy goes to protons 
before it reaches the proton gyroscale. However, when 
Pi ^ 1, heating of protons is marginal and most of the 
turbulence energy will cascade down further, crossing the 
proton gyroscale (see Quataert 1999; Quataert & Gruzinov 
1999). 

The Alfven waves will be converted into EMHD waves 
(whistlers) below the proton gyroradius scale. When 
EMHD turbulence is anisotropic (i.e. k\[ <^ k±), heating 
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of protons by EMHD turbulence will be marginal because 
"protons sample a rapidly varying electromagnetic field in 
the course of a Larmor orbit (Quataert & Gruzinov 1999)." 
In this paper, we confirmed that anisotropy of EMHD tur- 

1 /3 

bulence (fc|| cx /c/ ) is strong. Therefore, in this case, the 
remaining energy (that is, the energy that has survived the 
coUisionless damping before the proton gyroradius scale) 
cannot heat protons. Instead, it will heat electrons when 
it reaches the electron gyroradius scale. The bottom line 
is that the strong anisotropy of EMHD turbulence will 
make heating of protons by EMHD turbulence rather dif- 
ficult in coUisionless astrophysical plasmas with Pi 1 
(see Quataert & Gruzinov 1999). 

9.7. Other implications 

Other applications of the EMHD model include coUi- 
sionless magnetic reconnection in laboratory and space 
plasmas (Bulanov, Pegoraro, & Sakharov 1992; Biskamp, 
Schwarz, & Drake 1995; Avinash et al. 1998). Turbulence 
on microscale may be a source of anomalous resistivity, 
which can stabilize X-point reconnection. How impor- 
tant this type of reconnection is is a subject of debates. 
For instance, a model of magnetic reconnection in Lazar- 
ian & Vishniac (1999) appeals to the magnetic field weak 
stochasticity, rather than microphysical plasma effects to 
explain fast reconnection. Simulations by Kowal et al. 
(2009) successfully tested the Lazarian & Vishniac (1999) 
model, which may mean that in many astrophysically im- 
portant cases the reconnection is fast irrespectively of the 
plasma properties. The anomalous effects, however, may 
be important for the initiation of the reconnection when 
the original level of turbulence in the system is low. In ad- 
dition, in a partially ionized gas where the field wandering, 
which is the key element of Lazarian & Vishniac (1999) 
model, is partially suppressed, the anomalous plasma ef- 
fects can be important for reconnection (Lazarian, Vish- 
niac & Cho 2004). Incidentally, in the latter paper it is 
predicted that MHD turbulence is not killed by neutral 
friction in the partially ionized gas, but it gets resurrected 
at the scales at which neutrals and ions decouple. The 
resurrected cascade involves only ions and electrons, not 
neutrals, and may proceed as electron MHD cascade below 
the proton gyroscale. 

In addition, properties of EMHD turbulence are impor- 
tant for understanding of physics of neutron star crusts 
(Gumming et al. 2004; Harding & Lai 2004) and acceler- 
ation of particles in Solar flares (Liu, Petrosian, & Mason 
2006). While in most cases we view the EMHD cascade 
as the continuation of the MHD cascade below the proton 
gyroradius, for the crust of a neutral star the EMHD tur- 
bulence can be present on much larger scales. The main 



assumption of the EMHD is that one can ignore the mo- 
tions of protons. This is definitely the case of the neutron 
star's solid crust. 

10. SUMMARY 

We have found the following results. 

1. Electron MHD (EMHD) and Electron Reduced 
MHD (ERMHD) show identical scaling relations 
(both spectra and anisotropies). 

2. High resolution EMHD simulation confirms k~'^/^ 
spectrum obtained by earlier studies (Biskamp, 
Schwarz, & Drake 1996; Biskamp et al. 1999; Ng et 
al. 2003; GL04). The spectrum of electric field is 
consistent with k^^/^ spectrum obtained by earlier 
studies (see Schekochihin et al. 2009; Howes et al. 
2008a; Dmitruk & Matthaeus 2006). 

3. Our detailed study of anisotropy using different 

1/3 

techniques supports fey oc k_l the EMHD scaling 
obtained by CL04. 

4. Decaying EMHD turbulence and driven EMHD tur- 
bulence show the same scaling. 

5. When we use three or larger number of points to de- 
fine the structure functions, the scaling exponents 
of high-order structure functions follow a scaling 
similar to that of incompressible hydrodynamic tur- 
bulence. 

6. Bispectrum of ERMHD turbulence, reflecting the 
coupling of different scales in the cascade, looks very 
different from that of standard MHD one. 

7. The probability distribution functions (PDFs) of 
the increment of the magnetic field strength in 
EMHD and MHD cases are well described by the 

Tsallis distribution. The general trend of the PDFs 
is consistent with observations of the solar wind. 
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APPENDIX 

A. DERIVATION OF ERMHD EQUATIONS 

Schekochihin et al. (2009) first derived ERMHD from kinetic RMHD equations. They also gave derivation of ERMHD 
equations from the generalized Ohm's law. Therefore the starting point of ERMHD may be also the generalized Ohm's 
law. Derivation of ERMHD is almost identical to the EMHD case. 

Here we briefly summarize the derivation of the ERMHD equation from the generalized Ohm's law. Detailed derivation 
can be found in Appendix of Schekochihin et al. (2009). The magnetic induction equation reads 

— =Vx(v,xB), (Al) 



Table 1 
Runs 



Run 


Resolution 


Bo 


b at t=0 


at t=Ot 


fc|| at t=Ot 


Comments 




E512D 


512^ 


1 


1.21 


2 < fe < 5 


2 < fc < 5 


decaying 






E256D 


256=* 


1 


1.21 


2 < fc < 5 


2 < fc < 5 


decaying 






E256F 


256^ 


1 


1.21 


k ~ 2.5 


k ~ 2.5 


forced 






E256D-EL 


768x256^ 


1 


0.071 


4.5/^2 < fcx < 4.5v^ 


1/3 < /C|| < 1 


decaying 






ER256D1-EL 


768x256^ 


1 


0.071 


4.5/^2 < fcj_ < 4.5^2 


1/3 < /c|| < 1 


decaying. 


I = 


1 


ER256D8-EL 


768x256^ 


1 


0.07lt 


4.5/^2 < < 4.5^2 


1/3 < fc|| < 1 


decaying. 


1 = 


1/8 


MHD512F 


512^ 


0.8 





k ~ 2.5 


k ~ 2.5 


forced 






MHD256F 


256^ 


1 





k ~ 2.5 


k ~ 2.5 


forced 







^In this case, the value of b at t=0 is listed. 

tin cases of driven turbulence, the average driving-scale wavcnumber is listed. 
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Cho & Lazarian 



where 



Ve = V, - ^ = V, - —^V X B (A2) 

en„ 4,Tren„ 



(see Appendix C of Schekochihin et al. 2009). Here we ignored magnetic dissipation. Note that the usual EMHD equations 
are also based on these equations. Substituting Eq. (|A2p into Eq. (jAip . we get 

— ==Vx [v, xB--^(VxB)xB], (A3) 
ot Aneue 

= • - X X B) X B] - £ (vi-) X [(V X B) X B], (A4) 

« -BV • V, - — X [(V X B) X B], (A5) 

= -BV • V, - — ^— V X [B • VB], (A6) 

which is the same as the EMHD equation except the BV • term on the right. Here we assume « 0, but V • = 
V • Ve ^ 0. We ignore the last term of Eq. (|A4[) because it is ^ O(e^). According to the RMHD ordering (see Section 2 
of Schekochihin et al. 2009), we have 

— ^ 1^ ~ 7^ = e « 1, (A7) 

Tie Bq k± 



where Sn^ is the fluctuating density. 
For 6|| and bj^, we have 



where we use Vx « V_lX and BV • = BqV • + O(e^). Here e ~ k\\/k±. 
From the continuity equation, we can rewrite the BqV • (= BqV • Vg) as 



id \ Six 

BoV • V, = -Bo — + Vie • — (AlO) 

\Ot J He 

^ -B„(|+v,..V.)^^^^^^ (All) 
^^"+7^7T^^^^x(-^^x^ll^) (A12) 



/3i(l + Z/r) dt /3,(1 + Z/t 

2 9b 



m + Z/r) dt ' ^^^^^ 

where we use b|| — 5||Z, v^e x 6||Z = O(e^) « 0, and b\\/Bo — — (/3i/2)(l + Z/T){5ne/ne), which is equivalent to the 
pressure balance. Here, Z = qi/e (e = \qe\) is the ion-to-electron charge ratio, r — Ti/Te is the temperature ratio, is 
the mean electron number density, Sric, is the fluctuationing electron number density, Pi = SirnikBTi/ Bq is the ion plasma 
beta. Therefore, Eq. (IA9|1 becomes 

^ " -2 + m + Z/T)4nenJ^ ^ ' ^^^^^ ^^^^^^ 

In this paper, we use notations different from those in Schekochihin et al. (2009). The notations in this paper and in 
Schekochihin et al. (2009) are related by 

b ^ SB, (A15) 

(A16) 

bi 5BjL (A17) 
n-e ^ noe- (A18) 

B. TWO-POINT SECOND-ORDER STRUCTURE FUNCTION AND SPECTRUM 

We can obtain the 2-point second-order structure function of a variable A from EA{k): 

5F2(r) = < |A(x + r)-A(x)|2 > (Bl) 

= J d^xj d^k u(fc) e*-^(e* '' - 1) J d^k' u(fc') e'^' ''{e'^' '' - 1) (B2) 



EMHD Turbulence 



= J d^k u(fc) u*(fc) [2 - e"^ " - e-'^-"] 
^2Tr J dk k^\u{k)\'^ J de sin 61 [2 - e 



kr COS i 



- ikr cos 6 



TT / dk fc2|u(/c)|^ 



sin kr 



kr 



cx87rr" 



d{kr) (kr) 



knr 



1 



sin kr 
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(B3) 
(B4) 
(B5) 
(B6) 



where we assume 

/c2|u(fc)|2 oc fc"" (B7) 

for ko < k < oo. In the hmit of kr <C 1, the integrand is proportional to ~ (fcr)^^™. Therefore, if m < 3, we can rewrite 
Eq. dH]) as 

d{kr) {kr)^"" 



SF2{r) w STrr'^ 



1 



fcr 



If 771 > 3, however, the integral in Eq. (jB6|) is roughly proportional to (fcor')(/cor) 
In summary, we have 

r™-! ifm<3 



(fco?") and we have 



SF2{r) cx 



if m > 3. 



(B8) 



(B9) 



(BIO) 



In case of Kolmogorov, m = 5/3 < 3 and we have SF2{r) oc r"^ ^ — r^/"^. In case of EMHD, m for perpendicular direction 
(to=7/3) is still smaller than 3 and we have SF2{r) oc r™~^ = r'*/'^. However, that for parallel direction is expected to be 
larger than 3 if turbulence is anisotropic. Therefore, the 2-point second-order structure function for parallel direction is 
not suitable for revealing the true scaling exponent and we will have SF2{r) cx for parallel direction. 

C. MULTI-POINT SECOND-ORDER STRUCTURE FUNCTION AND SPECTRUM 

Falcon et al. (2007) and Lazarian & Pogosyan (2008) used the 3-point second-order structure function that can work 
with a steeper E{k): 

SF2{r) = < |A(x + r) - 2A(x) + A(x - r)!^ > 



= J d^xj d^k u{k) e*''-^(e* "" -2 + e"* '') J d^k' u(fc') e*'-^(e*' "" - 2 
= J d^k u(fc) u* (fc) [e*-'' - 2 + e"^''-'"] [e'*'-'" - 2 + e"*'''] * 
= 271 J dk k^\u{k)\^ J de sme[e 



r COS y I ^ — 2zkrcos0 -|- 6 ^g ikrcos^ — zkrcosf 



TT dk k^\u{k)f' 



sin 2kr 



2kr 



cx onr 



d{kr) {kr) 



knr 



3- 4 

4- 4 



-I- e' 
sin kr 



kr 
sin kr 
kr 



sin 2kr 
2kr 



(CI) 
(C2) 

(C3) 

(C4) 

(C5) 

(C6) 



where we assume k \\i{k)\ cx for ko < k < oo. In the limit of kr <C 1, the integrand is proportional to ^ (fcr) 
Therefore, the integral is roughly a constant (i.e. independent of r), if m < 5. If to > 5, however, the integral is roughly 
proportional to (fcor)(fco ?')'*"'" = {kor)^~"\ Therefore, we have 



SF2{r) cx 



„m — 1 



if TO < 5 

(fcor)^~™ cx r"' if TO > 5. 



Similarly, we can construct a 4-point second-order structure function: 

SF2{r) = < \A{yL + 3r) - 3A(x + r) + 3^(x - r) - A(x - 3r)p >, or 

SF2{r) = < |A(x + ^r) - 3A(x + ^r) + 3A(x - ^r) - A(x - ^r)|2 >, 



which scales as 



SF2{r) cx 



r™-i if TO < 7 

r'"-i(A:o?-)^"" cx if to > 7. 



We can also construct a 5-point second-order structure function: 

SF2{r) =< |A(x + 2r) - 4^(x + r) + 6A(x) - 4A(x - r) - ^(x - 2r)p >, 

which scales as 



SF2ir) cx 



^(fcoO^-™ cx if TO > 9. 



(C7) 

(C8) 
(C9) 

(CIO) 

(Cll) 
(C12) 



